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A POSTERIORI ERROR ESTIMATES WITH POINT 
SOURCES IN FRACTIONAL SOBOLEV SPACES 

FERNANDO D. GASPOZ, PEDRO MORIN, AND ANDREAS VEESER 


Abstract. We consider Poisson’s equation with a finite number of weighted 
Dirac masses as a source term, together with its discretization by means of 
conforming finite elements. For the error in fractional Sobolev spaces, we pro¬ 
pose residual-type a posteriori estimators with a specifically tailored oscillation 
and show that, on two-dimensional polygonal domains, they are reliable and 
locally efficient. In numerical tests, their use in an adaptive algorithm leads 
to optimal error decay rates. 


1. Introduction 


We consider the problem 

N 

(1.1) — Alt = ctjSxj in n, u = 0 on dfl, 

i=i 

where C is a polygonal domain with Lipschitz boundary dfl and, for any 
j = l,2,...,iV, each ajS^^ is a point source given by aj G M and the Dirac 
measure 6xj at xj € D. Such point sources are a useful idealization in modeling and 
appear in various applications: for instance, in modeling the effluent discharge in 
aquatic media [4] , in reaction-diffusion problems taking place in domains of different 
dimension [ 10 ], and in modeling the electric field generated by point charges [ 2 ]. 

Since D C K^, point sources induce singularities of the type log] • —Xj\. In 
particular, they do not belong to and so the solution u of ( 1 . 1 ) is not in 

Hence, the variational formulation within iJ^(D) cannot be used for (1.1) 
and the usual approach to a posteriori error estimation in the energy norm is not 
directly applicable. 

Although u ^ i7^(D), it can be approximated by finite element methods arising 
from the variational formulation within e.g., consider the finite element 

space V 7 - of continuous functions that are piecewise polynomial up to degree I over 
a triangulation T of D. Then the Galerkin approximation U € V 7 - given by 



is well-defined thanks to the continuity of the functions in V 7 -. 

In view of u ^ the error of the approximation U has to measured in a 

norm weaker than the H^-norm. For quasi-uniform meshes, Babuska [5] and Scott 
[18] derived a priori estimates for the error in the i/®-norm, s G [0,1). Exploiting 
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that the singularity of point sources is known, Eriksson [11] proves a priori estimates 
for the L^- and lE^’^-norms that show the advantage of suitably graded meshes. 

Several a posteriori error estimates, which may be used to direct an adaptive 
algorithm, are also available. Araya et. al. [3, 4] derived a posteriori estimates for 
L^’-norm, p S ( 1 , oo), and the lE^’^'-norm, p S {pq, 2 ), where po € [ 1 ; 2 ) depends on 
n. More recently, Agnelli et. al. [1] obtained a posteriori estimates for the weighted 
Sobolev norms introduced by D’Angelo [9]. 


In this article, we analyze a posteriori estimators for the error between u and 
its Galerkin approximation U in the i7^“®-norm, 0 < 0 < i. It is based upon 
the variational formulation within of Necas [15]. The error indicators are 

given by 


riT,e = 


II A17||^2(7') + 111^^^111^2(37') 


Ter, 


where Ht = stands for the local meshsize and |V?7] denotes the jump of the 

normal derivative across interelement edges. Our main result is the reliability and 
the local and global efficiency of these a posteriori error estimators. More precisely: 

Main Result. For any 0 < 9 < the error ofllG Yj- is bounded from above and 
below in the following manner: 

tgT 

1 

I?T.e < Ol I|u — 17 || 771 _ 9 (^^) {TgT), ( ^ < Cl\\u — U\\jji-e(^Qy 

T&r 


The quantity is an oscillation-type term (see Section 3.2) and the constants 
Cjj, Cl depend on 9, D, the polynomial degree i, and the minimum angle of the 
underlying triangulation. Moreover, ojt is the patch of all the elements sharing a 
side with T. 


The following comments are in order: 

• The indicators T £ 'T, coincide with the standard residual ones for 

Laplace’s equation and so do not depend on the point sources in (1.1). 

• The term (g, which depends on the point sources, vanishes under certain con¬ 
ditions, e.g., if the point sources for every star of the triangulation have the 
same sign; see Remark 3.2. Noteworthy, if this is not already given on the 
initial grid, it can be always met after a finite number of refinements. The 
term (g is thus a non-standard oscillation term. 

• In the many cases in which (g vanishes, the error of U is encapsulated only 
with the help of the approximate solution, without invoking data. 

The rest of the article is organized as follows: Section 2 reviews fractional Sobolev 
spaces, the variational formulation of ( 1 . 1 ) within , its dual problem, and its 
finite element discretization. In Section 3, we give a complete definition of the 
oscillation term f^g and prove the presented main result. Finally, in Section 4, we 
test the indicators rjT.e, T gT, in an adaptive algorithm. 

Throughout this article, a < b will denote a < Cb with some constant whose 
dependence will be stated whenever it is not clear from the context. We will use 
a ~ 6 to denote a < 5 and b < a. 
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2. Continuous, dual, and discrete problems 

In this section, we review fractional order Sobolev spaces, as well as the well- 
posedness of (1.1) and its dual problem in such spaces. This will be instrumental 
for building up the proofs of the a posteriori error bounds. Moreover, we recall the 
finite element discretization of (1.1) and associated notation. 


2.1. Fractional order Sobolev spaces. Let G be a bounded open set of d G 
N, with a Lipschitz boundary. We use the following notation for the (semi)norms 
of the usual (Hilbertian) Sobolev spaces H^{n) of integer order n G Nq: 


Iloilo, G ■“ \^\o,G 



u\u 




where, for any multi-index /3 = {Pi, 132, ■ ■ ■, 13d) € Nq, its length is defined by 
1/31 '■= Pi + P 2 + ■■■+ Pd and denotes the weak /3-derivative of p, with the 
convention D^^’---'^'>p = p. If s > 0 is not an integer, we write s = n + t with n G Nq 


and 0 < t < 1 and define the norm of iJ'^(n) by ||<(>||g := \\p\'^ 


n,G 


l</>l 


s,G’ 


with 


I<(>L.g — 


E 

|/3|=n 


\DPp{x) - Df^p{y)\''^ 


G JG 


\x-y\ 


d-\-2t 


dx dy 


The space Hq{G) is the completion of C^{G) under the norm || • ||s_g. 

The following lemma summarizes some basic properties of the fractional Sobolev 
spaces iL®(n), which can be found, e.g., in [14]. 


Lemma 2.1 (Fractional Sobolev spaces). Let G be a bounded open subset 
d G N, with a Lipschitz boundary. We have: 

(i) (Shift) If \P\ < s, then p G H^{G) entails D^^p G iJ®-l^l(G). 

(ii) (Trace) If s > there exists a constant G such that 

(2.1) ||</>||o.aG := ||</'||L^(aG) < C'||</.|U,G, for alipGG^{G), 

and thus the trace operator, defined in C°°{G) astrp = p\do can be extended 
to be a continuous operator from (G) into L^ {dG). In particular, we have 
trp = 0 for all p G Hq{G) and s > |. 

(iii) (Sobolev embedding) If s > ^ then II‘^{G) is continuously embedded into 

G^’*(G), with k = [s— — 1 and s— | = k + t. More precisely, if s—^ = k + t 

with fc G No and 0 < t < 1, then for any function p G H^{G) there exists a 
function p G G*(G) such that p = p a.e. in G and 


max II LoofGl + max sup 
l/3|<fc ^ ^ m = kx,yeG 

x^y 


\D^Pix)-DPp{y)\ 
\x - y\* 


<C\\P\\s,G- 


The constants G appearing in (ii) and (iii) depend only on the set G, the dimension 
d, the smoothness order s, but are otherwise independent of p. 
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We shall need several inequalities involving fractional Sobolev seminorms on 
domains (bounded, connected and open sets). The first one is the counterpart of 
the classical Poincare inequality for functions in 

H^G) := |(/)e : J = 0 

We shall use it to derive further inequalities, as well as for the proof of the lower 
a posteriori error bounds. 


Lemma 2.2 (Fractional Poincare inequality). // 0 < s < 1 and G is a bounded 
domain ofW^, there is a constant Cp depending on s and G such that 

mlG<Cp\^\lG. forallgbGH^G). 


Proof. The simple proof of Faermann [12, Lemma 3.4] (which readily generalizes 
to d-dimensional domains) shows 


( 2 . 2 ) 


^ 1 diam(G)^+2" 

- 2 |G| 


where diam(G) and |G| stand for the diameter and the Lebesgue measure of G, re¬ 
spectively. The dependence on s can be identified more precisely, see, e.g., Bourgain 
et. al. [8], but (2.2) suffices for our purposes. □ 


We need also the following generalization of the Friedrichs inequality. It is related 
to the well-posedness of (1.1) and the definiteness of the error notion considered in 
the a posteriori analysis below. 

Lemma 2.3 (Fractional Friedrichs inequality). Let s > 1/2 and G C &e a 
bounded domain with Lipschitz boundary dG. We have 

(2.3) UWo.G < CF\f>UG, for all ef G H^{G), 

where the constant Gp depends on d, s and G. 


Proof. We distinguish different cases for s. 

Q If s G N, the claimed inequality is the Friedrichs inequality for Sobolev spaces 
of integer order; see, e.g., [7, Ch. II, 1.7]. 

[ 2 ] Let I < s < 1 and fix ^ G Hq{G). For any constant c G K, we write 

where ji9Gj denotes the {d — l)-dimensional Hausdorff measure in Thus, the 
Cauchy-Schwarz inequality on dG and the trace theorem (2.1) imply 

1/2 

||/>||o.G < !!/> - cJlo.G + ^ 11^ “ cJlo.SG < !!</ - cJlo.G + C\\(j) - cJls^G 

with G depending on G. We choose c = [G]”^ (f and obtain the claimed inequal¬ 
ity in this case with the help of Lemma 2.2. 

[ 3 ] Let k < s < k + 1 with /c G N and assume, without loss of generality, that 

(j) G G((°(G). Observe that, for any multi-index /3 with j/3j = fc, we have = 0 

and, therefore. Lemma 2.2 yields 

\\D^f}\\o,G < C\D<^^\s-k,G, whence [(/[^.g < C\(I)\s,g- 

The claimed inequality then follows from Step [j~]. □ 
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For any s > 0, we define H ®(G') as the topological dual space of Hq{G). In 
other words: is the set of linear functionals 'ip : Hq{G) —)■ M satisfying 

|(V^,</))| <C||<^|U.G, forall</>ei7J(G), 


for some constant G independent of (p. Hereafter, the symbol (•, •) indicates the 
result of the application of a functional to a function in its domain of definition. 

The spaces H~^{G) with s > ^ are of particular interest for us. It is convenient 
to equip them with norms that have simple scaling properties, like seminorms. 
Lemma 2.3 ensures that, for such s, the seminorm | • |s_g is a norm in Hq(G), 
equivalent to || • ||s,g- We thus can define a suitable norm for H~^{G) with s > 5 
by duality: 


11, II W’,0) 

\m-s,G ■= sup J- - = sup J- -, 

<peH‘{G) ms,G 'PeC^iG) ms,G 


for iP G H-^{G). 


We then have also 

(2.4) I<^L,g= sup for G iLo(G')- 

For the sake of brevity, we will write iL®, H^, || • ||s, to denote 7 Jq(H), || • 

respectively. 

Finally, we need a Bramble-Hilbert-type inequality; it will be useful in deriving 
the a posteriori upper error bound. For this purpose, it is sufficient to consider 
triangular domains. Moreover, in view of the following lemma, it is sufficient to 
derive such, and other, inequalities for the reference triangle T in given by the 
vertices ( 0 , 0 ), ( 1 , 0 ), and ( 0 , 1 ). 


Lemma 2.4 (Scaling properties of Sobolev norms). Let T C LI be a triangle, set 
hx := and let F(x) = Ax + b be an affine bijection such that F{T) = T. 

(i) Assume that the functions (p : T and (p : T satisfy <p = <po F. Then, 
for any s > 0, we have (p G F[^{T) if and only if (p G H’^{T) and 

(ii) Assume that the distributions ip : G^{T) —>■ M and ip : G^{T) —>■ M satisfy 
{ip,poF) = I det(H)|“^('0, (/)) for all p G G§°{T). Then, for any s > ^, we 
have Ip G H~^{T) if and only if ip G H~^(T) and 



The hidden constants depend only on s and the minimum angle of T. 


Proof. Up to constants depending on the minimum angle of T, we have 

|det(H)|~hp, ||H|| ~ hr, and ||H“^|| ~ 

The case s G Nq is well known, see, for instance, [7, Ch. II, 6 . 6 ]. Next, consider 
0 < s < 1. Given p G H^{T), we derive 


\€,t = 




\p{x) - Pi'yffi 


< 


det(H )|2 f f \p{x)-p('y)\^ 


) det(H)p dxd'jj 


p-l||-(2+2s) |^_^|2+2, 


■ dx dy hip ^101 


2 

s,f’ 
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where we have used that \A{x — y)\ > ||v4“^||“^|a; — y\. The opposite inequality 
follows in a similar manner using \A{x — y)\ < ||^|||i — y\ and (i) is thus verified 
also for 0 < s < 1. 

For the case fc < s < fc + 1 with fc S N, we observe 


^ \DP4>{.x)-Df^<l>{y)\^\\A-^t Y, \DPm-D^^{x)\ 

\p\=k \p\=k 


and then proceed for each term on the right as in the case 0 < s < 1. 
It remains to verify (ii). Given ip G H~^{T) with s > we obtain 


-s,T = 


sup T7|- 

<t>(iC^(T) W\s,T 


sup 

06C’o“(f) 


\ det{A)\{ip,^) 


,1-s 


l</>L 




iv^l- 


s,T 


with the help of (i). 


□ 


Lemma 2.5 (A fractional Bramble-Hilbert inequality). Let 1 < s < 2. There is 
constant Cbh depending only on s such that, if (p G H‘^(T) vanishes at the vertices 
ofT, then 

Proof. We let Ii denote the Lagrange interpolation operator onto the polynomials 
7^^ of degree at most 1. Given any polynomial p G of degree at most 1, we may 
write 

(P = (P-ii(P= {(p - p) + Ii{(p - p). 

Since f < Cmaxj. \ip\ for some constant C, Lemma 2.1 (hi) yields 

II 0 II 1 J ^ C\\<p-p\\^ f, 

where C additionally depends on s > 1. We choose p G such that Jfp = J.p(p 
and Jfdip = JfdiCp for i = 1,2. Thus, we can conclude with the help of the 
classical Poincare inequality and its counterpart Lemma 2.2. □ 

2.2. Continuous and dual problem. We assume that 11 C is a two-dimen¬ 
sional polygonal, but not necessarily convex domain with Lipschitz boundary 511. 
For any 9 > 0, Lemma 2.1 (iii) then implies that is continuously embedded 

in C°(ll) and, therefore, the right-hand side of Poisson’s equation in (1.1) belongs 
to the space H~^~^. 

Writing s = 1 — 0, we are thus led to consider the Dirichlet problem 

(2.5) — Au = / in n, m = 0 on 511, 

where / G and s < 1. Expecting that the solution is in H^, we additionally 

require s > ^ to provide a meaning to the boundary condition in (2.5) by means 
of Lemma 2.1 (ii). Furthermore, since we shall invoke duality arguments, it will be 
useful to consider also the range 1 < s < |. In any case, the differential operator 
—A should be understood in the distributional sense: 

{—Au,ip):={u,—Aip), V(^ G C(j"(H). 

Integration by parts shows that, if u G Cq°{LI), then 

(2.6) {v , — Aif) = B[v, ip] := / Vv-Vipdx. 

Ja 
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Necas [15] extended ;B to a continuous bilinear form on Hq x Hq ^ by means of the 
following lemma. 

Lemma 2.6 (Weak derivative and fractional spaces). Let 0 < 6 < 1/2 and G be a 
bounded Lipschitz domain o/K‘^. Then there exists a constant C such that 

f ^9dx<C\\f\U_s,G\\g\kG: forallf,g&C^{G). 

JG OXi 

Noteworthy, the case 0=1/2 has to be excluded in view of Grisvard [13, Propo¬ 
sition 1.4.4. 8 ]. Moreover, [15] verified that the extension of B on x satisfies 
an inf-sup condition and thus obtained the following theorem. 

Theorem 2.7 (Weak formulation in Hq). Let | < s < |. For any f G 
there exists a unique weak solution u G Hq of the Dirichlet problem (2.5); 

B[u,q,] = {f,^}, V^GCo°°(fl). 

It satisfies the a priori estimate 

lhL<c|i/L_2, 

where the constant G depends only on s and fl. 

Theorem 2.7 may be used to establish the well-posedness of Problem (1.1) in 
suitable fractional Sobolev spaces. 


Corollary 2.8 (Well-posedness for point sources). Problem (1.1) has a unique 
solution u which, for every 0 > D, satisfies u G H/^~ and 

N 

ld[u,(l)] = '^aj(j){xj), y(j) G 
i=i 


Let us fix 0 < 0 < i. In view of the extension of the bilinear form B to 


^ X 77^+® and Corollary 2.8, we refer to 

given / G H~^~^, find u G such that 

B[u,c/] = (/,</.)_!_,,1+,, V</.G i7i+^(fl) 

as the direct (primal) problem and to 

given g G , find w G such that 

B[ip, w] = {ip, , Wif G 

as the adjoint (dual) problem. 


(2.7) 


( 2 . 8 ) 


2.3. Finite element discretization. Let T be a conforming (edge-to-edge) tri¬ 
angulation of n. We refer to the minimum angle appearing in T as shape coefficient 
ap. Moreover, we denote by £ the set of all its edges and by V the set of all its 
vertices. The star around a vertex z G V is given by 

w. := U T. 

TeT-.zeT 

Given 7 G N, we let denote the space of polynomials of total degree < 1. 
Moreover, let V 7 - be the finite element space of continuous piecewise polynomials 
that vanish at the boundary, i.e. 

Vr = = {P G C(n) : P = 0 on dfl; P|r G iP^VT G T}. 
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The set of the standard nodes (the locations of the degrees of freedom) of Yj- is 
indicated by Afi. We thus have V n C A/^, with equality for £ = 1. 

The approximate solution U of (1.1) is defined as the Galerkin solution in V 7 -: 

N 

(2.9) [/ G Vr : B[U, V]=J2 ajV{xj) VV € Yr- 

Notice that, although Problem (1.1) is associated with a non-symmetric weak for¬ 
mulation, U is computed by solving the usual symmetric, positive definite linear 
system. 


3. A POSTERIORI ERROR ANALYSIS 

In this section, we derive a posteriori bounds for the error \u — U\i- 0 , where 
0 < 6 * < 1/2 and | • |i_e is a norm thanks to Lemma 2.3. Before embarking on their 
derivation, we define estimator and data oscillation. 

3.1. Estimator and data oscillation. We start by defining the error estimator. 
To this end, we denote by hr = |T |2 the local meshsize and let E G £ be an edge. 
If E is an interelement edge, we write E = Ti D T 2 with Ti, T 2 G T and define the 
jump |VC/]|£; of the flux by 

lYUj\E ■= VC/i • -b Vl/2 • 

where U^, denote the restrictions of 17 to Ti, T 2 , respectively, and n^, are 
the outer normals of Ti, T^- If E is a boundary edge, we have E C d£l and set 
|V[/1|£; := 0. With these notations, we define estimator and indicators by 

(3.1) ife = ^ 7 t,6i ?7T,e = ll^t^l!L 2 (T) + III^^lllL 2 (aT) ■ 

tgT 

Notice that rjTfi depends only on the approximate solution U and the local meshsize 
hr and, thus, is independent of the point sources in Problem (1.1). 

The oscillation is tailored to the specific class of source term in Problem (1.1). 
It is defined starwise and depends on the interplay of the boundary, the nodes and 
the supports of the point sources. We write 

77/ = Aft U d£L 

for short. For each vertex z G V, we consider only point sources whose supports 
are not nodes and collect them according to their sign: 

At ■■= {J : Xj G uJz\ Aft and aj > 0}, := {j : Xj G uJz\ Aft and aj < 0}. 

Given j G A+, we set 

cr/" := min < dist(a;j,7i/)® + max dist(a:i,7i/)® , max \xj — Xif >, 

I iGAJ leA- j 

if A~ 7 b and atj = 0 otherwise. Similarly, given j G A~ 

cr J ■ := min < dist(a:j,77/)® -b max dist(a:i,7i/)® , max \xj — xtf >, 

I *gA+ ieA+ j 
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if A'^ ^ 0, and tT“j = 0 otherwise. Moreover, let Az denote the piecewise affine 
function that is 1 at z and 0 in all other vertices of the star Wz. The oscillation 
indicator associated to z is 


^eiz) 


dist{xj,J^e)^\aj\Xz{xj) if z e dfl, 

min { Ejga+ II 


with the convention ^0 = 0. The global oscillation term is then 


(3.2) 


^0 '■= 



2 


Let us conclude this section with two remarks, which highlight useful properties of 
this oscillation. 

Remark 3.1 (Scaling of oscillation). The error estimator r]g involves, as scaling 
factors, suitable powers of the local meshsizes hx, T G T- For example, the L^- 
norm of element residual At/ is scaled by h}fX^. In order to derive a bound for 
with similar scaling factors, we observe that, for any vertex z of a triangle T G T, 
we have 


where the hidden constant depends on the shape coefficient ax- Since the triangles 
of the patch ojt = : T' n T 7 / 0} have similar areas, this readily leads to 


(3.3) 


^0 ^ 



,|T(r) 


Xj Gcjt 


with 


T(r) = 


if T n = 0 and aiUj > 0, Vxi, Xj G ujt, 
otherwise. 


The scaling factor of the oscillation is thus given by hfp. In the case of quasi¬ 
uniform refinement, the global counterpart of this scaling factor corresponds to the 
decay rate of the error \u — U\i-g under consideration; see [18]. The bound in 
(3.3) may be also used as a triangle-indexed alternative for (3.2), which however 
overestimates whenever the distances between point sources are much smaller that 
the local meshsize. 


Remark 3.2 (Oscillation and refinement). The local oscillation ^g{z), z G V, 
vanishes in many cases and, in any event, asymptotically. To see this, we observe 
that ^g{z) 7 / 0 implies, according to the location of z, one of the following conditions: 

• If z S dn is a boundary vertex, there is a point source located in \ Afe- 

• If z S 0 is a interior vertex, there are point sources with different sign located 
in Wz \ Mi- 

Remarkably, if one of these condition is verified, a finite number of refinements step 
ensures that its negation is met and this remains so for further refinements. Hence, 
= 0 can always be reached after a finite number of suitable refinement steps. 
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3.2. Upper Bound. The estimator (3.1) and the oscillation (3.2) provide an upper 
bound for the error in . 


Theorem 3.3 (Upper bound). Let u he the solution of Problem (1.1) and U its 
approximation associated with the triangulation T■ There exists a constant Cu, 
depending on LI, 9 G (0, |) and the shape coefficient aj- of'T, such that 

\u - C^li_e < Cuive + Ce)- 

It is worth observing that this upper bound will simplify under adaptive refine¬ 
ment. 

Remark 3.4 (Asymptotic form of upper bound). Point sources generate singu¬ 
larities which are centered at their supports. Since an adaptive algorithm will refine 
around these places, Remark 3.2 suggests that, after a finite number of adaptive 
refinements, the upper bound of Theorem 3.3 becomes 

\u - U\^_g < Cuve- 

This expectation is in line with our numerical experiments in §4. Interestingly, the 
asymptotic form is independent of the point sources in Problem (1.1). 

We now prove Theorem 3.3, postponing some technical estimates about the 
interpolation error to Lemma 3.5 below. 


Proof of Theorem 3.3. We split the proof in several steps. 

nn We start by relating the error with a suitable norm of the residual. To this end, 
we use a duality argument relying on 


\u-U\^.e 


sup 

ggiZ-l + O 


{u-U, g) 

ll<?ILi+. 


which is a special case of (2.4). Let g G H and denote by w the solution to 
the dual problem (2.8). Theorem 2.7 and the original symmetry of B in (2.6) yield 
that w G with Q blLi+e- Since 


{u — U , g) = B[u — U, w] 


we obtain 


N 

{R , w) := ajw{xj) — B\U, re], 


(3.4) \u-U\^_g<ClM\-i-e- 
It thus remains to bound the residual norm ||i?||_i_e. 

0 We rewrite (i? , w) for a fixed w G by means of the partition of unity 

~ ^ exploiting (i?, Xf) = 0 for all interior vertices z G V n fl. 

To this end, we let I-pw denote the Lagrange interpolant of w onto and set 
w := w — I-pw. Moreover, if x € V H H is an interior vertex, we let G M to be 
chosen later, while, if z G V n 11 is boundary vertex, we set Cz := 0. Then 

(3.5) {R, w) = {R, w - Irw) = '^{R, {w - Cz)Xz) 

zGV 


with the local contributions 


(3.6) {R, {w - Cz)Xz) = ^ aj[w{xj) - Cz]Xzixj) - B[U,{w - Cz)Xz] 
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and Az = U . 

Fix any z € V. Assuming that = w{x^) for some G uJz, we bound the 
second term in (3.6) as follows: 

\B[U,{w-Cz)X.]\<( Y, 4+®||A;7||o^^+ ^ 4^®||lVC/l||o^^ 

Eee^ 

where He indicates the length of an edge E and Tz ■= {T G T '■ z G T} and 
£z '■= {E G £ : z G E} stand for the triangles and interior edges of the star 
respectively. To this end, we integrate by parts on each T G Tz, use A 2 < 1 and 
obtain 




|S[C7,(iZi-c.)A.]|< ^ \\AU\\,^^\\w-Cz\\,^^+ Y IIIVC/1IIo,eII^-c.|Io,e- 
TgT, Eg£. 

Adopting standard arguments to the setting at hand, see Lemma 3.5 below, yields 

\\w - CzWq e ^ and \\w - Cz\\q e ^ \Mi+e,u>,, 

where the hidden constants depend on 6 and the shape coefficient < 77 -. Inserting this 
in the preceding inequality, we arrive at the claimed bound for \B[U, (w — Cz)A^]|. 
0 Let z GV . For given or appropriately chosen c^, we derive the following bound 
for the sum over Az in (3.6): 


(3.7) 


Y Olj[w{Xj) 


s] A2 (^Xj ) 


< 


?e(z) \w\ 


i+e,uj^ 


The possible choices of Cz depend on the location of z as well as the point sources 
located in uiz- 

Case 1: z G £l and Az = 0. Then ~ <^z\Xz{xj) = 0, irrespective of 

the choice of Cz- In particular, (3.7) is verified and we may take Cz = 0. 

Case 2: z G and Az ^ 0. Then we have A~ 0 or A+ ^ 0. If the latter occurs, 
then J2kGA+ OLkXz{xk) > 0 and we can consider 


which implies 


(3.9) Y = Y 

Fix j G A~ for a moment. On the one hand, the definition of cj" and Lemma 3.5 
below yield 


(3.10) \w{xj) - c+l < Y 

iGA+ 


< max 
i^At 


Xj-X^\ \w\i+g^ui. 


On the other hand, since w{v) = 0 for all v G Afe, we have that, for any choice of 
Ui GUt.iG {0}U A+, 


\w{xj) - c+l < \w{xj) - w{vo)\ + Y 

i&At 
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which again by Lemma 3.5 implies that 

(3.11) — c+l < ( dist(a:i,+ maxdist(a:i,./\/f)® ) |w|i+e . 

V i&At ) 

Combining the bounds (3.10) and (3.11) with the definition of a~ gives 
\w{xj) - c+l < a~\w\i+e,oj,, for all j G A~, 
and thus, upon recalling (3.9), 

(3.12) ^ aj[w(a;j)-c+]A;,(a;j) < ^ \aj\a~X^{xj)\w\i+g,^^. 

Moreover, if we have A~ 0, we can consider 

e: = y: w.th ,3r ,= ^ > 0 

Notice that the definition of c~ is the one of c+, if we replace A~ by A+. Conse¬ 
quently, we can argue as before and obtain here 

(3.13) ^ aj[ic(xj) - c^]A^(a;j) < ^ \o:j\cr^Xzixj)\w\j^^g^^^. 

Taking the minimum of the two bounds (3.12) and (3.13) verifies (3.7) in this case. 
Case 3: z G dXl. Here we have Cz = 0. Using Lemma 3.5 and w{x) = 0 for all 
X G Afe, we derive 

- Cz]Az(xj) = ajw{xj)Xz{xj) 


jeA, 


jeA, 


< 


Y \aj\dist(xj,JCtf ) k| i+e,u 

yieA, 


which verifies (3.7) also in this case. 

Notice that all choices of Cz in Cases 1-3 satisfy min,^^ w < Cz < max;^^ w. Since w 
is continuous and compact, we can always choose x^ G uJz such that w{x^) = Cz 
and therefore apply Step 3. 

0 Combining Steps 3 and 4 and using He hr whenever E C T, we derive 


{R, {w- Cz)Xz) 


< 


Y + ^e{z) 


.TeE 


\w\i+e,uj. 


We insert this inequality into (3.5), use the Cauchy-Schwarz inequality for sums 
twice, observe that the cardinality of Tz is bounded in terms of and arrive at 

\ 1/2 


{R, w) < I 

VzGV 


li+e.wj 


Subsequently, X^zgv l^li+e^ — (3-4) finish the proof. 

We turn to the postponed estimates about the interpolation error. 


□ 
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Lemma 3.5 (Interpolation error). Let w € with 0 < 9 < ^ and consider 

w := w-I-tw, where I'fW denotes the Lagrange interpolant ofw intoYij-. Moreover, 
let Wz, z € V, he any star ofT- Given any x,y S ujz, we have 

\w{x) - w{y)\ < |x - yf\w\i+e,u:,, 
and, if Cz = w{x^) for some G Wz, then 

\\w — CzW^ rp < hjf \w\i+e^u)^, W'^ ~ ^z\\o,E 

for any triangle T G T and any edge E G £ containing z. The hidden constants 
depend only on LI, 9, I and a^, while He stands for the length of E. 


Proof. We start by deriving the first inequality, where the domain is the refer¬ 
ence triangle T instead of some star Thanks to Lemma 2.1 (iii), the bound 
f ^ Riaxj, |w|, and Lemma 2.5, we have, for x,y GT, 


|w(x) - w(y)| 
\x - y\^ 


< 




where the hidden constant depends only on 9 and £. In view of Lemma 2.4 (i) 
and its proof, both sides scale in the same way under affine transformations of the 
domain. Hence, for any triangle T G E, we obtain 

(3-14) \w{x) - w{y)\ < |x - yf\w\i+e,T- 


If X and y are two arbitrary points in we connect them with a polygonal path, 
made of straight segments in each element of lUz and having total length < \x — y\. 
The existence and construction of such a path is presented in Lemma 3.4 of [17], 
where the involved constant depends on uj- and the Lipschitz constant associated 
with dLl. Applying (3.14) segmentwise, we obtain the first claimed inequality. 
Integrating it, we readily deduce the other ones. □ 


3.3. Lower Bounds. In this section, we assess the sharpness of the upper bound 
in Theorem 3.3, dealing with the two parts rje and fe separately. 

Let us start with the oscillation fg defined in (3.2). We first recall that, typically, 
oscillation terms are not shown to be bounded by the error, but are, formally, of 
higher order. Here, we encounter similar properties for fg. Indeed, Remark 3.1 
suggests that, under global uniform refinement, fg decreases at least with the order 
of the error jit — Lf\i+g. Moreover, Remark 3.2 suggests that f^g even vanishes after 
a finite number of appropriate refinements in a reasonable adaptive algorithm. In 
such cases, fg is then of arbitrarily higher order. 

Our main result about the sharpness of rjg from (3.1), or of the asymptotic form 
of the upper bound in Remark 3.4, is as follows. 

Theorem 3.6 (Lower bounds). Let u be the solution of Problem (1-1), U its ap¬ 
proximation associated with the triangulation E and 0 < 0 < 1/2. For any triangle 
T G E, we have the local bound 

VT ,0 < Cl \u- , 

where ujt is the patch of all triangles ofE sharing a side with T. Furthermore, we 
have also the global bound 

Ve ^ Cl\u — U\i-g 
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Both constants Cl, Cl depend only on 9, the polynomial degree I, the shape coef¬ 
ficient aj- and the number N of point sources. 

Remark 3.7 (Asymptotic independence on point sources). The dependence on 
N is actually through the maximum number Nj- of point sources supported in one 
element of T. After a finite number of suitable refinement steps, every triangle will 
contain at most the support of one point source. For the same reason as in Remark 
3.4, one expects that these refinement steps are actually quickly accomplished by 
a (reasonable) adaptive algorithm. We therefore may say that the constants Cl 
and Cl are asymptotically independent of N. Combining this with the asymptotic 
upper bound in Remark 3.4, we see that, asymptotically, the error \u — U\i-e is 
encapsulated with a posteriori quantities that are independent on the point sources 
in Problem (1.1). 

The proof of Theorem 3.6 uses the constructive approach of Verfiirth [19]. In 
order to adapt it to our setting with fractional Sobolev space at hand, we need 
the following preparations concerning the local continuity of B and suitable test 
functions with local support. These test functions will be products of polynomials 
and cut-off functions. For the latter, we shall use the following type: given a ball 
of radius r with midpoint z € K^, set 

where ,W=(“»(ld^) ^ 

1^0 otherwise. 

Lemma 3.8 (Cut-off within triangles). Let k gNq and T be a triangle. Moreover, 
let B be the ball with maximal radius in the reference triangle T and F’ : —>■ 

be an affine bijection with F{T) = T. Then the cut-off function riT := rj^ o F~^ 
satisfies, for all V S V^{T), 

[v^<fv^r,T and {VvtMo.t < \\V\\^ ^ . 

JT JT 

The hidden constants depend only on 6, k, and the minimal angle ofT. 


9b{x) = p 


X — z 


Proof. In view of the transformation rule and Lemma 2.4, the claim is equivalent 
to 


< 


i/2 


Vb 


and \y{Vpg)\gf< 


V 


0,T 


for all V G V^{T). This statement in turn follows from the equivalence of norms 
on the finite-dimensional spaces 7^^(r) and 7^^(T’)/K. □ 


Lemma 3.9 (Cut-off across edges). Let fc £ No and A = Ti C T 2 be the common 
edge of two triangles Ti and T 2 . For i = 1,2, denote by T) the reference triangles 
with vertices (0,0), (1, 0), (0, (—1)*) and indicate by B 12 the ball with maximal radius 
in the reference patch Ti U T 2 . Moreover, let F : IR 2 ^ piecewise affine 

bijeetion with F{Ti) = T^, i = 1,2. Then the cut-off function pE '■= 9 b ° 
satisfies, for all V G V^{E) and i = 1,2, 

[ V^< [ PE and \y(VpE)\e,T,<hp-<^\\VpE\\,E<h~E'^~'\\y\\oE^ 

J E J E 

where Le denotes the length of E and V is a suitable extension of V. The hidden 
constants depend only on 9, k, and the minimal angle of T. 
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Proof. In view of the transformation rule and Lemma 2.4, the claim is equivalent 
to the following statement associated with the reference edge given by the vertices 
(0,0), (1,0): for all V € V^{E), we have 


< 




Vb 


and 




VVb 


< 


0,Ti 


V 


0,E ’ 


where the extension of 14 is given by 14(aii, 0 : 2 ) = t4(a;i). Again, these inequalities 
follow from the equivalence of norms on the finite-dimensional spaces V^iTi) and 


Lemma 3.10 (Local continuity of B). Given 0 < 9 < 1/2 and any triangle T, we 
have 

J Vv ■ Vipdx < \v\i-g^T\'^ip\e,T 

for all V G and g) € such that suppt/J = suppiyr or supp77£: with 

It, Ve from Lemmas 3.8 and 3.9. The meaning of the left-hand side is given by 
continuous extension and the hidden constant depends on the minimal angle ofT. 


Proof. In view of Lemma 2.4, both sides of the claimed inequality scale in the same 
manner under affine transformations. We therefore can assume that T = T, where 
T is one of the reference triangles T, Ti, T 2 . Correspondingly, we write B for B or 
Bi 2 . We set c = and apply Lemma 2.6 to obtain 

/ '^v-\7ipdx= / \7{v — c) ■ dx < \\v — c\\.^_g f\\\7ip\\g f. 

JfJf ' ’ 

Using Lemma 2.2 and replacing dG by the 2-dimensional set T\B in step 2 of the 
proof of Lemma 2.3, we conclude with 

Ik — c||i_6i,T ^ and ||V:/5||6i,r ^ D 


After these preparations, we are ready to prove the claimed lower bounds. 


Proof of Theorem 3.6. [//I We shall use test functions, whose support does not con¬ 
tain point sources. To construct them, we shall exploit Lemmas 3.8 and 3.9 for the 
following sub-triangles. Let Nj- be the maximum number of Dirac masses supported 
in a triangle T gT and set 

(3.15) M ■.= 2{Nr + l). 

We divide each edge if of T into M equal sub-edges and denote by the set of 
these sub-edges. Moreover, for any triangle T G T, we join the endpoints of the 
sub-edges by lines parallel to the edges of T and so divide T into equivalent 
sub-triangles. We indicate with Sff this set of sub-triangles; see also Figure 1. The 
choice (3.15) of M ensures that 

(a) for each T G T, we can choose a sub-triangle T* G that does not contain 
any point source appearing in (1.1), 

(b) for each edge B of the triangulation T, we can choose a sub-segment E* G 

such that no point source appearing in (1.1) is supported in the union 
of the two sub-triangles adjacent to E*. 
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Figure 1. Sub-edges, sub-triangles and support of test functions. 

Scaling arguments, similar to those in the proofs of Lemmas 3.8 and 3.9, yield: 
if fc € No, we have 

(3.16) l|l/|lo,T<l|l"llo.T. and \\V\\,^e < \\V\\o,e^ 

for all V G 7^^(T) or F € V^{E), where the hidden constants depend on M, k, and 
the minimal angle in T but not on the choices of T* and E*. 

0 Let us now prove a lower bound of the local error in terms of any given element 
residual I1 ^N||q j,, T G T. To this end, we only need to consider £ > 2 and 

observe that then AC7 G V{TY~‘^. Using tpr = rjT-, we derive 

\\AU\\Ie<\\^U\\Ie^< [ AU^t = - f VU-VipT= f V{u-U)-VipT 
Jt* jt* Jt* 

|u — j,, |V:/?T|e,T* Y: II^HIIo.t 

with the help of the first inequality in (3.16), Lemma 3.8, integration by parts, the 
choice of T*, Lemma 3.10 and Ht- < hT < MHt*. Multiplying by h}jf^ ||A[/||Qi^ 
and squaring we arrive at 

(3.17) hl+^<^\\AU\\l^<\u-U\l_,^^. 

0 Next, we provide a lower bound for the local error in terms of any given jump 
residual |![VC/]||q E edge of T. Let Ti, T 2 be the two triangles of T sharing the 
edge E, and let , U^, n^, denote the restrictions of U and the outer normals 
of Ti, T 2 , respectively. Then |V[/] = G V^~^{E) and denote 

by |V17] its extension from Lemma 3.9. Using ipE = [VU] tje*, we derive 

WUjg^E= Y. f 

i=l,2 

\7{U — u) ■ VipE + AU pe 

< 1^ ~ + II AU||q J.. ||v5e||o,T* 

1 = 1,2 

< ^ h-^-"|n-t7|i_,,,J|[V[/l||o,E + 4l|At/||o^^J|[Vt/l||o.E 
1 = 1,2 

with the help of the second inequality in (3.16), Lemma 3.9, integration by parts, 
the choice of E*, Lemma 3.10 and He* < hE < MhE-. After multiplying by 


= Y, Si'U ■ \/pE + AU Pe = 


i=l,2 ■ 


i=l,2 ' 


livc/lll 


0,E 


< 


livuii' 


0,E* 


< 
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ll[VC/]|!o_^, we obtain 

(3.18) 4^®lllVC/l||o_^< ^ 

2 ^1,2 

m The claimed lower bound in terms of rix,T G T, follows by combining the squares 
of (3.17) and (3.18) for the involved triangles and interelement edges; recall that 
fiE ^ It-t whenever E is an edge of T and that for boundary edges E C dfl, we 
have set |V17]|£; = 0. 

The global lower bound is a direct consequence of local one: sum the square of 
all local ones and take into account that the cardinality of {T' G T \ ujt' G) T} is 
bounded in terms of the shape coefficient of T. □ 

4. Numerical Results 

In this section, we numerically test the a posteriori error estimators of §3. To 
this end, we use it in the adaptive solution of two examples of Problem (1.1) and 
analyze resulting properties of the adaptive algorithm. 

The adaptive algorithm, which was implemented within the finite element tool¬ 
box ALBERTA [16], has the following structure. Given 6 and a conforming initial 
triangulation To of 17, it iterates the main steps 

(4.1) Solve — >• Estimate —Mark — >• Refine. 

The step Solve consists in solving the discrete system (2.9) for the current trian¬ 
gulation T and linear elements. The step Estimate then computes the a posteriori 
error estimator (3.1) and the step Mark selects triangles for refinement by means 
of the maximum strategy: T gT is marked whenever > 0.5 max^/gT-g. In 
the step Refine, these marked triangles are bisected twice so that each of their 
edges is halved. In doing so, further triangles are bisected in order to maintain the 
conformity of the next triangulation. 

Example 4.1 (Fundamental solution). Consider Problem (1.1) with data such 
that 

=-ii-logkl, x e 17 := (-1,1)^ 

ZTT 

is the exact solution, together with the parameter values 9 — f), 0.125, 0.250, 0.375, 
and 0.5 for the adaptive algorithm. Notice that, for 0 = 0, the error estimator 
formally corresponds to the infinite error \u — U\i^q. Moreover, 9 = 0.5 is not 
covered by the analysis given in §3, but the convexity of (—1,1)^ suggests that the 
equivalence of 771/2 and \u — C/|i/ 2 ,n is still given. 

For given 0 > 0, the exact solution u formally has almost 0 derivatives (in 
L^) more than required in the error \u — U\i-e^Q. Thus, with increasing 0, quasi¬ 
uniform meshes ensure increasing error decay, suggesting that the grading of meshes 
generated with rjo decreases with increasing 0. Figure 2 confirms this expectation, 
as well as the corresponding meshes for the intermediate values 0 = 0.125, 0.375, 
which are not shown. If 0 is small, the mesh grading is very strong: for 0 = 0, the 
triangles at the origin of meshes with about 1000 degrees of freedom (DOFs) have 
areas smaller than 10“^®. In the case of 0 = 0.125, this happens for meshes with 
about 5000 DOFs. 
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Figure 2. Mesh grading and error norm for the fundamental solu¬ 
tion. The triangulations after 20 iterations of (4.1) for 9 = 0.000, 0.250, 
0.500 (from left to right) illustrate that the mesh grading decreases with 
increasing 6, which corresponds to weakening the error norm. 


A next step in our numerical testing of the estimator rjg could be to study 
the decay rate of the estimated error \U — u|i_6/,n- This will be done for the 
second, more involved example. Here we shall instead study decay rates of two 
error notions, for which the estimator rjg is not originally designed. The first error 
notion is \u — U\i^qo with = fl \ H( 0 ; |) = {ix,y) S : |a;| + |j/| > j}. 
Since u G the maximum decay rate for it with linear finite elements is 

^DOFs”^^^, reached for example by uniform refinement. The second error notion 
is the L^-error \u — 17|o,n. Here we have D'^u G LP{^1) for any p G (0,1) and thus 
u is an element of the Besov space Bp{LP{^)). Consequently, the maximum decay 
rate with linear finite elements is #DOFs“^, reached for example by thresholding 
[ 6 , Theorem 5.1]. Figure 3 suggests that, except for 0 = 0, the meshes generated 
with r]g provide the maximum decay rate for both errors. We however notice an 




Figure 3. Convergenee histories of the -error off the singularity 
(left) and the L^(Jl)-error (right) for the fundamental solution. In both 
cases, we plot error versus #:DOFs in logarithmic scales. The plots for 
0 = 0, 0.125 and 0.25 prematurely end because we encounter triangles 
whose areas are below 10~^®. 

advantage for greater values of 9, where the initial stagnation of the error decay is 
shorter. This stagnation expresses the difference between the observed error notion 
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and the estimator, which puts more importance to the singularity at the origin. For 
0 = 0, the fact that u ^ H^{n) appears to be reflected in an infinite stagnation. 

Example 4.2 (Point sources and reentrant corner). Consider the non-convex L- 
shaped domain ft = (—1,1)^ \ ([0,1) x (—1,0]) and the boundary value problem 

— Am = 5 ( 0 . 33 , 0 . 66 ) + 5(_o.251.-0.85) + 5(_o.25,-0.87) in ^7 
M = 0 on d^l. 

The exact solution is not known to us; see Figure 4 for an approximate solution 
and corresponding triangulation. Notice that boundary condition and right-hand 



Figure 4. Approximate solution and triangulation for Example 4.2 
after 25 iterations with 9 = 0.375. 

side are compatible at the reentrant corner. Nevertheless, the refinement at the 
reentrant corner indicates some non-compatibility of the error. The three point 
sources entail that the exact solution has three logarithmic-type singularities, two of 
them being very close. As for Example 4.1 and independently of data compatibility 
at the reentrant corner, we thus have D^u G L^(n) for all p G (0,1). For 9 G 
( 0 , 1 ], the maximum decay rate for the approximation with linear finite elements 
in I • |i_e,n is therefore again reached for example by thresholding 

[ 6 , Theorem 5.1]. Figure 5 indicates that these maximum decay rates are also 
obtained with the triangulation generated with the help of the estimator rjg ; machine 
precision prevents a possibly better confirmation in the cases 9 = 0.125 and 0.25. 
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DOFs 


e 

slope 

i-\-e 

2 

0.125 

-0.484 

-0.562 

0.250 

-0.612 

-0.625 

0.375 

-0.684 

-0.688 

0.500 

-0.748 

-0.750 
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